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ABSTRACT 

The Boltzmann model for the random generation of "decomposable" combinatorial 
structures is a set of techniques that allows for efficient random sampUng algorithms 
for a large class of families of discrete objects. The usual requirement of sampling 
uniformly from the set of objects of a given size is somehow relaxed, though unifor- 
mity iimong objects of each size is still ensured. Generating functions, rather than the 
enumeration sequences they are based on, are the crucial ingredient. 

We give a brief description of the general theory, as well as a number of newer 
developments. 

1 INTRODUCTION 

1.1 Random Generation of Combinatorial Structures: a 
Quick Overview 

Random generation is often used as a tool for exploration ("what do large objects of 
this type look like?"), or to provide large datasets for software and algorithm testing. 

Usually, one defines a combinatorial class as a finite or countable family of 
discrete "objects", equipped with a "size" function from to the natural numbers, 
such that for each natural number n, the subset 'rfn of all objects with size n is finite. 
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A uniform random generator for the class is then a randomized algorithm that takes 
as input an integer n, and outputs an element randomly selected uniformly from '^n- 
The efficiency of the generator is typically measured in terms of its expected time and 
space complexities, expressed as a function of the size n. 

Random generation methods tend to fall into one of a small number of classes. 
Ad hoc methods rely on precise combinatorial properties of the considered class; a 
fine example is provided by Remy's algorithm [14] for the random generation of plane 
binary trees. Markov chain methods rely on the simulation of a Markov chain whose 
states are the objects one wants to sample from, and that converges to the uniform (or 
another suitably chosen) distribution; careful analysis of the convergence speed makes 
it possible to run the chain for a fixed number of steps and obtain an almost-uniform 
generator, or more sophisticated techniques such as Coupling from the Past [13] can be 
used for exact sampling from the stationary distribution. 

Decomposition methods apply for classes where objects are, informally, "made up" 
of smaller objects, be they from the same or another class that is itself decomposable. 
A prime example is that of plane binary trees, that is, rooted trees where each internal 
node has exactly two children, one of which is distinguished as the left child while the 
other is the right child. In this case, a plane binary tree is either made up of a single 
root-leaf, or of a root and left and right subtrees, both of which can be any plane binary 
trees. 

The first systematic example of a decomposition method for random generation is 
the so-called recursive method [11], where decompositions are used to obtain recur- 
rences satisfied by the counting sequences. These counting sequences are then used, 
together with the decomposition rules themselves, to guide the random generation al- 
gorithm. In the above example of plane binary trees, the counting sequence (the well- 
known Catalan numbers) is used to determine the probability p„^k that a uniform ran- 
dom tree with n intemal nodes has a left subtree with k internal nodes, and to sample 
K from this distribution; then the random generator is recursively called with sizes K 
and n — 1 — A" to obtain the left and right subtrees (which are independent conditioned 
on their respective sizes); the resulting tree is then uniform. 

1.2 The Boltzmann Method 

The Boltzmann method, as introduced in [7, 8], is another decomposition-based method 
that can be applied to roughly the same combinatorial classes as the recursive method. 
We will give a precise description of the method in the next section. For this introduc- 
tion, we will simply describe the crucial ingredients. 

Where the recursive method uses counting sequences for the random generation 
algorithms, and generating functions are mostly a tool to compute the counting se- 
quences, the Boltzmann method uses the generating functions, viewed as analytic func- 
tions of a real variable (and, in practice, the values of the generating functions) in the 
random generation algorithms - thus reducing the need for precomputation to a small 
number of real constants. The idea is to "relax" the requirement for a uniform sam- 
pler (which outputs a uniform random structure among those of the target size n) into 
allowing structures of all sizes, while keeping uniformity among all objects of each 
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individual size. By choosing the "right" distribution on sizes, independence among 
substructures is introduced, which results in very simple and efficient algorithms. 

1.3 Outline of the Paper 

Sections 2-4 make up the bulk of what can be termed the "Boltzmann method". Sec- 
tion 2 gives a general description of the Boltzmann method. Section 3 lists a number 
of constructions which can be used to define specifications of combinatorial classes for 
which Boltzmann samplers can be automatically compiled from the specification. Sec- 
tion 4 sums up various results on the complexities of the random sampling algorithms. 
Section 5 describes how Boltzmann samplers can be used to get closer to the classical 
model of uniform, fixed- size random generation. 

Section 6 deals with the question of how one can effectively obtain the real con- 
stants used in Boltzmann samplers, and Section 7 describes results that draw on the 
principles of the method without exactly fitting in it. 

The list of references does not attempt to give a complete list of articles using the 
ideas exposed here. The interested reader will find more examples in the bibliography 
of [2], even though it is not limited to references about the Boltzmann method. 

All proofs and most technical details have been purposefully omitted; the algo- 
rithms in Section 3 have been included mostly to demonstrate their simplicity. We 
have made the choice of not detailing any of the many examples that could be given; 
the interested reader will find many such examples, including pictures of large random 
structures, in the original papers. 

2 THE BOLTZMANN METHOD: GENERAL 
DESCRIPTION 

The Boltzmann method can be used in two flavors, the ordinary (or unlabelled) and 

exponential (labelled) variants. 

Throughout the paper, we use the word structure (or '^i' -structure) to indicate an el- 
ement of a combinatorial class 'lo. No particular assumption is ever made on the nature 
of such structures, though classical examples tend to come from discrete mathematics 
or theoretical computer science: words over some finite alphabet, sequences, various 
flavors of trees, etc. 

2.1 Combinatorial Classes and Products 

Let 'To be some combinatorial class. The size of an object c e wiU be noted |c|. For 
any integer n, let c„ denote the number of objects in ^ with size n. The ordinary (resp. 
exponential) generating function for is 

C(z)=£c„z«, resp.C(z)=£c„^; 

it is always assumed that the considered generating function has positive radius of 
convergence p, i.e. that limey" < oo (resp., lim(c„/n!)^/" < oo). 
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For < X < p, the (normal, resp. exponential) Boltzmann distribution over for 
parameter x is the probability distribution defined, for any c e by 

fJc) = , resp. YJc) = 

' C(x)' ^ ' n!C(x) 

These distributions give positive probability to all objects in the class, with the 
property that two objects with the same size have the same probability. They are, 
of course, not the only probability distributions with this property; their interest lies 
mostly in their relationship with two common constructions in the combinatorial world: 
the normal and labelled products. 

Given two classes £/ and their normal product is just their Cartesian product 
= ^ X ^, with size defined additively by | (a, b) \ = \a\ + \b\ for (a, b)e£/x^. 

To define the labelled product, one has to assume that structures are made up of 
both an unlabelled structure c and a "labelling", a permutation on a set whose is that of 
c. Think of a structure c as being formed of |c| basic "atoms", each of which receives 
a distinct label from [[1, |c|]] = {i G Z : 1 < / < |c|}. By a slight abuse of notation, we 
identify these atoms with the integers 1 to |c|, so that the labellings are just permutations 
C7 e ■ For s^ch unlabelled structure c, the set of admissible permutations may be a 
strict subset of J^^^.^ . 

Then, the labelled product of two labelled structures (a, (7a) and {b,Gb) is defined 
as {a,b), with admissible labelings obtained by taking all partitions of [[1, |a| + \b\]] 
into two parts A and B of respective sizes \a\ and \b\, and, for each partition, taking 
the one permutation a € '^\a\+\b\ where all entries in A are in the same respective 
order as that of Oa, and all entries in 5 are in the same respective order as that of o^, 
(that is, if (x,y) € A^, then (7o(x) < ab(y) iff a{x) < C7(y), and similarly for B). As 
a result, (a,(Jo) and {b,ai,) have ('"'q'*') different structures in their labelled product. 
The labelled product of two classes is defined as the set of all labelled products of 
structures in the two original classes, with size again defined additively. 

The first important property is as follows: if two unlabelled (resp., labelled) classes 
£/ and 3i have generating functions A (z) and B(z) (resp., exponential generating func- 
tions A(z) and B{z)), then their normal (resp. labelled) product has generating function 
C{z) =A{z)B{z) (resp., C{z) = A(z)B(z)). 

An immediate, and most useful, consequence, valid under both models, is this: 
= £/ xM, then taking the product of two independent si^ -structure and SS-structure, 
each following the Boltzmann distribution with parameter x, results in a ^-structure 
under the Boltzmann distribution with parameter x. For labelled structures, it is implied 
that one selects a uniform random set to define the permutation in the product. 

This "independence under substructures" property, in turn, has interesting practi- 
cal consequences, in that it makes it possible to describe, for a number of classical 
combinatorial constructions, systematic ways to produce efficient algorithms to sam- 
ple from the Boltzmann distribution for classes that are entirely described (possibly in 
a recursive way) with them. This is the topic of Section 3. 
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2.2 The Boltzmann Method 

By a Boltzmann sampler for a combinatorial class "if, we mean a randomized algorithm 
Vrf that takes as input a real parameter x and outputs a random "^-structure under 
the Boltzmann distribution with parameter x. Our overall goal is to create efficient 
Boltzmann samplers for as many combinatorial classes as possible. 

Given a description of a combinatorial class 'rf from which we would Uke to obtain 
"large random structures" (of size n, ideally), the Boltzmann method can be summa- 
rized as follows: 

1. Find out if our class can be specified (up to a reasonably simple size-preserving 
bijection) from the constructions in Section 3. If not, try to extend the expressive 
power of the method by adding new constructions. If this does not work, maybe 
the Boltzmann method is not the best choice after all. 

2. Use the techniques of the Purple Book [10] to locate the "dominant" singularities 
(those of smallest modulus, which govern the asymptotics of the counting coeffi- 
cients; since the generating functions have nonnegative coefficients, at least one 
such singularity lies on the positive real axis) of the generating functions in our 
specification, and possibly an estimate of the value we should give parameter x 
to give expected size n to "^-structures under the Boltzmann distribution. 

3. Compute approximations (to roughly 0(logn) digits) of the values at x of all 
involved generating functions, possibly using the combinatorial oracle of Sec- 
tion 6. 

4. Use the patterns in Section 3, together with our specification, to write a Boltz- 
mann sampling program. 

5. Optionally, add a rejection scheme to obtain samples with size in [(1 — e)n, (1 + 
e)n], or even of exact size n (more costly). 

Alternatively, step 2 can be replaced by experimentation using the other steps. 



2.3 Choosing the Parameter and Tuning for Size 

As said above, the Boltzmann model uses a real parameter that may be chosen arbi- 
trarily inside the radius of convergence of the generating function for the class under 
consideration. This parameter x governs the distribution of sizes of the random struc- 
tures. In typical applications, one would like to obtain large structures, that is, there is 
an ideal value n for sizes. We now briefly turn to the question of picking an appropriate 
value of X for a target n. 

For a given x, the probability that a Boltzmann-distributed ^ -structure will have 
size n is given by 

/ \ afjX^ 

multiplying by n, and sunaming over all values of n, we obtain the expected size 

/ , xA'lx) 
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where A! is just the derivative of the generating function A (which can be formally 
defined as A! (x) = na„x"^ ' ; this corresponds to the usual derivative inside the radius 
of convergence of A). 

Except in degenerate cases, Nj;/ is a strictly increasing and convex function on the 
interval [0,p), and the equation N{x) = n has at most one solution; if, as is often the 
case, A'(x) goes to infinity as x goes to p, the equation has a unique solution x„ for 
each integer n. Setting x to this value x„ in the Boltzmann samplers from the previous 
section results in a sampling algorithm that produces structures of expected size n. 
Interestingly enough, the same equation also describes the value of x that maximizes 
the probability /?„ that the output structure will have size n. 

When the generating function is known exactly, one can solve for the exact value 
of x„. In many situations of interest, the generating function is known only through an 
equation that it satisfies (this is typically the case when the class is defined recursively, 
as described in the next section). In this case, it is often possible to use the techniques 
of analytic combinatorics [10] to derive precise asymptotic information about the gen- 
erating function and obtain a precise estimate of x„. 

In some cases, an attractive alternative is to use the singularity x = p; although it 
often implies an infinite expected size, a simple adaptation of the sampling algorithms 
makes this a very viable choice. This will be described in more detail in Section 5. 

3 BASIC CONSTRUCTIONS AND BOLTZMANN 
SAMPLERS 

In this section, we describe a number of constructions that can be used to describe 
more complex classes from simpler ones, and, for each construction, the corresponding 
combination algorithm that allows one to build a Boltzmann sampler for the new class 
using Boltzmann samplers for the classes involved in the description. The constructions 
described here allow one to describe combinatorial classes that are close to those of the 
theory of combinatorial species [1]. 

When nothing is specified, these constructions apply to both labelled and unla- 
beled structures; in the labelled case, it is silently assumed that one performs a label 
redistribution as in the case of the labelled product. 

The initial constructions were described in [8]; later additions are credited individ- 
ually. 

In all cases, the construction is translated into an expression for the generating 
function of the new class in terms of the previous one; this in turns gives a simple 
construction for the sampling algorithm, where "substructures" are independent. Many 
algorithms can be expressed in the form "Draw integer k from discrete distribution 
li(x), then let / receive the concatenation of k independent calls to generator Ts^ (x)"; 
we abbreviate this as 

y^[rM(x) =^ r^(x)]. 

We also use samplers for a few standard distributions: Bernoulli with success prob- 
ability X (Bern(x)), geometric (with support N) with parameter x (Geom(x)), and Pois- 
son with rate x (Poiss(x)). A subscript condition on these samplers means a condition- 
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ing on the output, which can be achieved by rejecting outputs until the condition is 
met. 

3.1 Finite Sets 

Finite (typically small) sets do not require an elaborate theory to produce sampling 
algorithms, and are included to serve as elementary bricks for more complex construc- 
tions. 

One typically defines an "empty structure" class S', containing a single structure of 
size zero that we denote as 1, and an "atom" class jF, containing a single structure Z 
of size 1 . Occasionally, one may use a number of different atom "types", which will 
then be written Za,Zi„ and so on. 

3.2 Disjoint Union 

The most basic construction is that of disjoint union: if and ^ are disjoint classes, 
then their union = (with size inherited from the original class) is a new class, 

whose generating fimction is simply C{z) = A{z) +B{z). 

Algorithm r[^U^] 

if Bern(A(x) /(A(x) +B{x))) then 

Return Ts^ (x) 
else 

Return T^{x) 
end if 



3.3 Product 

The product construction, being of fundamental importance, has been described in the 

previous section. 

Algorithm T[s^ x ^] 
Return (r^/(x),r^(x)) 



3.4 Sequence 

If ^ is a class with no structures of size 0, ^ = Seq(j2/) is the set of sequences 
(Ai , . . . , A;t)j for arbitrary A: > 0, with A, e s^, and size defined additively by 

\{A,,...,Ak)\ = \A,\+--- + \Ak\. 

The generating function for ^ is the pseudo-inverse of that of 

and the corresponding Boltzmann sampler for is as follows: 
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Algorithm rSeq(^) 

7^ [Geom(A(x)) =^ r^(x)] 

Return y 



3.5 Cycle (labelled) 

If ^ is a class with no structures of size 0, = Cycle(^) is the set of cycles of 
stmctures, that is, sequences defined up to a circular permutation of the component 
^ -structures. 

Working with labelled structures means that each sequence of k structures has ex- 
actly k—1 other structures that correspond to the same cycle. Consequently, the gen- 
erating function is 

C(z)=I^ = -log(l-A(z)), 
k>i 

and the corresponding Boltzmann sampler uses the "logarithmic" distribution Hx{k) = 
k\\og(\-x)\ ^ 1)' provided by sampler Loga(): 



Algorithm rCycle(^) 

7^ [Loga(A(x)) =^ Ts^ix)] 

Return 7 



3.6 Set 

If ^ is a labelled class with no structures of size 0, ^ = Set(,G/) is the class of sets 
of ^ -structures, that is, (possibly empty) sequences up to an arbitrary permutation of 
component structures. 

The corresponding generating function is 

fc>0 

and the distribution for the number of components is the Poisson distribution: 



Algorithm r[Set(^)] 

7<(- [Poiss(A(x)) r^(x)] 

Return 7 



A set construction for unlabelled structures was introduced in [9]. It is based on the 
multiset construction described next, and requires more elaborate manipulations. 



3.7 Multiset (unlabelled) 

Boltzmann samplers for the Multiset construction were introduced in [9]. 
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If is an unlabelled class with no structures of size 0, ^€ = MSet(^) is the class 
of all multisets of -structures - sets with possible repetitions. 
The corresponding generating function is 



The Boltzmann sampler uses a MaxIndex(A,x) subroutine, which samples from the 
discrete distribution defined by 



The Boltzmann sampler itself is then as follows: 



Algorithm rMSet(^) 

k *r- MaxIndex(A,x) 
for j from 1 to — 1 do 

7 ^ 7 U [Poiss (A [x^ )/j) =^ j copies of Ts^ {x-i ) ] 
end for 

7^ 7U [Poiss>i(A(x*)/A:) ^copies of r^(x*)]. 
Return 7 



(It should be noted that k in the above algorithm is no f the maximum number of 
repetitions of a structure in the output, but only a lower bound: all calls to the sampler 
are independent, so that a structure may be output more than once and obtain larger 
multiplicity.) 

Note that, in contrast to the previous constructions, the value C(z) is expressed not 
in terms of the value A(z), but of the values of A for a whole geometric sequence of 
values. 

3.8 Recursive Constructions 

All of the above constructions can be used recursively, i.e. a class can be defined 
using one of these constructions on a class that is itself (ultimately) defined in terms of 
't^ itself. Some care must be taken to avoid circular definitions: recursive specifications 
define structures from smaller structures, possibly of the same type, but not from them- 
selves. Thus, one can define a class by ^ = ^ x Seq(^) (this defines plane trees: 
a plane tree is composed of a root having an ordered sequence (possibly empty) of 
children, each the root of a plane tree), but a specification such as ^ = x ^ would 
be invalid if class ^ contains structures of size (it would then attempt to create an 
infinite number of ^-structures of size 0). 

Subject to this "well-foundedness" condition, all the previous constructions can 
be used recursively - and indeed, in most applications of interest recursivity is used. 
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Whenever it is the case, the Bohzmann samplers derived from the previous subsec- 
tions become recursive algorithms, for which termination can only be guaranteed with 
probability 1 (and in finite expected time; see Section 4). 

An effective characterization of this "well-foundedness" condition is given in [12], 
for specifications that outright forbid structures of size 0. 



3.9 Ordered Structures and Differential Operators (labelled) 

The constructions described in this subsection appear in [15, 5]. 

The derivative a' of a labelled combinatorial structure a is obtained by replacing 
the atom in a having the largest label with a "hole" - this hole holds the place of an 
atom, but does not contribute to size and does not get a label. Thus, the derivative 
of a structure of size n has size « — 1. The derivative of a combinatorial class sz^ is, 
of course, the class ^ = s^' = {a! : a ^ £/} of all derivatives of ^-structures. The 
corresponding (exponential) generating functions are related by 

C{z) = A'iz), 

A{z) = ao+ [ C{z)dz, 
Jo 

where ao is the number of £/ -structures of size 0. 

Derivative classes can be used in recursive constructions, under suitable "well- 
foundedness" conditions (see [2, 5] for details), to define a class as the solution to 
a symbolic differential equation. In very rough terms, this corresponds to imposing 
order conditions on labels. A simple example is provided by the class ^ of decreasing 
(labelled) binary trees, that is, labelled binary plane trees where each node (atom) is 
required to have a larger label than each of its children: the equation reads 

% = 0. 

and should be understood as this: the largest label in a decreasing binary tree has to 

be at the root, so its derivative will either be a unique object of size zero or equivalent 
(after relabelling) to a pair of binary trees, each of which has to be decreasing. 

Bodini etal. [5] describe a generic Boltzmann sampler for a class defined by a first- 
order differential operator jz/' = ^(iF, .2/), provided one has a Boltzmann sampler for 
the class s^) (that is, ^ is defined in terms of other classical constructions, and 

the whole sampler will necessarily be recursive). Like the sampler for multisets, it 
requires a change of the parameter - here, a random change - for each recursive call. 

Given the generating function A and a parameter < xq < Pa, one defines a proba- 
bility density on the interval [0, 1] by 

xoA'(iao) 

if J7 is a random variable following this distribution, J/xo can be interpreted as the 
result of picking a random point (according to Lebesgue measure) in the domain < 
y < A(x) , < X < xo, and keeping the abscissa x. 

With this definition, the Boltzmann sampler is as follows: 
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Algorithm Ts/ from He/', ^' = ^(^,^) 
ifBem(A(0)/A(x)) then 

Return a uniform object from ^ 
else 

Draw U G [0, 1] following density 

Return object (Z,/) with atom Z having the largest label 
end if 



3.10 Multivariate Models 

So far, we have only considered generating functions with a single variable, which 
"counts" for the size of the structures; that is, each structure yin the class contributes a 
single term z'^' or z'^'/y! to the generating function. 

Given a combinatorial class one can define an arbitrary number of statistics s,- : 
^ ^ N (1 < J < A:), and the corresponding multivariate generating function C(z, ui,...,Uk) 
(first as a multivariate formal power series, then as an analytic function) by changing 
the contribution of each structure c S to z''^' ni<K*:M-'^'^' (for the rest of this subsec- 
tion, we assume unlabelled structures), and consequently, a Boltzmann distribution for 
any tuple {x,u\, . . . ,Uk) of positive real variables lying inside the convergence domain 
of the generating function. 

When the considered statistics are transmitted additively under the constructions 
described in this section (say, if there are several types of "atoms", and statistic Si 
counts the number of atoms of type 0, the Boltzmann samplers can be adapted to this 
generalized model. 

This area has not been explored as extensively as others, probably because a general 
theory would involve additional technical details. Bodini and Ponty [4] use it to sample 
from context-free languages with a nonuniform distribution where the frequency of 
letters is artificially skewed, with an application to "Tetris tesselations" (perfect tilings 
of a rectangular region with pentominoes) where each piece has the same frequency. 

4 ALGORITHM COMPLEXITIES 

Each of the individual algorithms in Section 3 has low overhead complexity, but they 
tend to make possibly unbounded numbers of calls to other algorithms. The general 
theorem below is a compilation of results from [8] and other papers that extend the 
expressive power of "specifiable" classes. 

Theorem 1 Let ^ denote a (labelled or unlabelled) class that can be entirely specified, 
in a possibly recursive way, with the constructions of Section 3, and let < x < p be 
any positive real lying inside the convergence domain for the generating function of 
^. Then, assuming an oracle that provides values of the relevant generating functions 
at real values, the algorithm Y^^ compiled from the specification by the patterns of 
Section 3 terminates with probability 1 and in finite expected time, outputs a random 
^-structure under the Boltzmann distribution with parameter x, and uses a nimiber of 
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real arithmetic operations that is linear in the size of the output. 



5 APPROXIMATE AND EXACT SIZE SAMPLERS 



The user of random generation algorithm is often used to exact random samplers (al- 
gorithms that take n as input, and output a uniform random element of the subclass 
'^S„), or perhaps approximate-size samplers (algorithms that take two integers n < N as 
input, and output a random element with size in [n,N] with equal probability for any 
two elements with the same size). Both types can be obtained by adding a rejection 
mechanism to Boltzmann samplers, at overall costs that depend on the distribution of 
sizes under the Boltzmann model. 

We denote IJ.i{x) for the expected size, Hiix) for the expected squared size, and 
a{x) for the standard deviation on size, aU as functions of parameter x; analytically. 



5.1 Approximate Size Samplers 

Assume we are given a target size n, some tolerance e > 0, and the value x„ of the pa- 
rameter that ensures that the expected size of structures is n. We obtain an approximate 
size sampler with acceptable sizes in / = I{n,e) = ((1 — e)n, (1 + e)n) by repeatedly 
using the Boltzmann sampler r^(x„), until a structure whose size lies in / appears. 

Since the cost of each call to the Boltzmann sampler is linear in the size of the 
output by Theorem 1, the expected cost of this approximate size sampler is asymptotic 
to n times the expected number of calls to the Boltzmann sampler. 

In favorable situations, described as "lumpy" distributions in [8] and characterized 
by a(x)/;Ui(x) — > as x — > p , this expected number of calls goes to 1 as « (the 
target size) goes to infinity - the probability that the first call will yield a structure 
whose size is in / is asymptotically 1. More precise information on the distribution 
of sizes produced by such an approximate- size sampler can be obtained through the 
asymptotics of the generating function. Typically, for "lumpy" distributions, this size 
is concentrated around the expected size, and shorter intervals of length o{n) could be 
used without altering the theoretical results. 

In less favorable situations, this success probability goes to a finite positive constant 
as n goes to infinity, so that the expected cost of the approximate size sampler is still 
asymptotically linear. In many cases, a size tolerance of, say, 5% around the target size, 
at the cost of a constant number of rejections, is quite acceptable. 

In some even less favorable situations, it may be necessary to change the class 
by using "pointing" (a combinatorial operation close to derivation, corresponding to 
distinguishing a single atom in the structure) a finite number of times before one gets 
to the situation described above. 
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5.2 Exact Size Samplers 

Setting e = 1 /n in the approximate size samplers of the previous subsection, results 
in exact size samplers. It should be noted, however, that in most cases, the success 
probability (the probability, with parameter x„, of obtaining a structure of size exactly 
n) is only of order 0(l/n), which results in an expected complexity 0(n^) for the exact 
size sampler. 

5.3 Singular Samplers 

Whenever the generating function is finite at its dominant singularity p, one can de- 
fine a Boltzmann distribution forx = p, and the Boltzmann samplers can be used with 
parameter p. This is typically (though not universally) true with recursive specifi- 
cations, and the most frequent case is for the generating function to have a "square 
root-type" singularity, i.e. C{z) has a singular expansion of the form C{z) = C{p) + 
a{l-z/py^^+o{{l -z/p) 1^) as z approaches p. 

In such cases, the expected size for the singular Boltzmann model is infinite. While 
this offers the best chances of success for the approximate size samplers, it also im- 
plies that the expected cost of the same approximate size sampler is infinite, which is 
unacceptable. 

However, on closer examination, this infinite expectation only comes from those 
(rare: the probability is ®{\/^Jn)) runs of the sampler where the output size is much 
larger than the target n. A simple modification of the Boltzmann samplers can thus 
avoid this higher cost, by keeping track of the number of atoms generated so far and 
aborting the Boltzmann sampler as soon as the total becomes larger than n; with this 
modification, expected costs for a square root singularity become 0(n) for approximate 
size with finite £, and 0(n^/'^) for exact size. 

6 GENERATING FUNCTION EVALUATION AND 
PRECISION 

One of the key points of the Boltzmann method, when compared to the recursive 
method, is that enumeration sequences are replaced by generating function evalua- 
tions. In many of the constructions of Section 3, each involved generating function 
needs to be evaluated for the same value of its variable; in a few of them, some have to 
be evaluated for a deterministic or random sequence of values. 

In some cases, the generating functions have closed forms and this evaluation does 
not lead to special complications, but recursive specifications lead to generating func- 
tions that are determined by equations, and the question of determining a good approx- 
imation of the required values becomes more troublesome. 

Pivoteau, Salvy and Soria [12] provide an efficient solution to this problem, at 
least for the basic constructions of sums, products, sequences, cycles and sets. The 
natural idea of iterating the equations provides only slow convergence; the preferred 
method is based on Newton iteration, which ensures quadratic convergence (asymptot- 
ically, distance to the exact solution is squared by each iteration). For a given value 
of the variable, the generating function equations typically have several real solutions. 



6 GENERATING FUNCTION EVALUATION AND PRECISION 



14 



only one of which corresponds to the generating function that is the "true" solution; 
an important result in the above-mentioned paper is that, thanks to the existence of a 
combinatorial equivalent to the Newton iteration, convergence to the "true" value is 
ensured. 

Also, note that, when applying Newton iteration for generating function evalua- 
tion, convergence is significantly faster the further x is from the dominant singularity 
p. While this is good for constructions such as multisets or ordered structures that 
require the use of sequences of values (for sequences of values of the variable that 
quickly decrease to zero), it is conversely bad news for classical applications that re- 
quire very large structures, since this means taking values of x that are very close to 
p. Nevertheless, the experimental results reported in [12], even for very complex spec- 
ifications, remain within reasonable bounds (the computation time for the oracles of 
specifications implying 500 equations, for values of the parameter leading to expected 
structure sizes in the tens of thousands, are of the order of a minute). 

Another question that arises naturally is that of the influence of approximations on 
the final distribution of random samples. Even assuming a "perfect" source of ran- 
domness for the simulations, a small error in a generating function value that is used 
repeatedly by a sampling algorithm might result in a significantly biased distribution 
of the final samples. 

To give an example, suppose a specification involves a disjoint union = si/U^, 
and each of ^ and ^ is defined recursively using ^ - this is not an artificial assump- 
tion. Assume that, for the value x of the parameter, A(x) is slightly overestimated, and 
B{x) is sUghtly underestimated. In this case, each time the sampler for is used, it 
will have a tendency to switch to ^ more often than it ideally should; this will result 
in a distribution that is biased in favor of "^-structures that often use ^ -structures as 
components. 

One possible solution [6] to estimate, and possibly correct, this bias, is to design 
Boltzmann samplers that not only output a random structure, but also, for each real- 
valued constant A ~ A(x) used in the sampling, a "safety interval" [A",A+], with A^ < 
A < A+, with a precise meaning of "if the sampler had been run with any value in 
[A^ , A+] instead of A, the result of the whole computation would have been the same". 
This is done by studying the small number of discrete distributions one really needs to 
sample from. 

From such "safety interval" samplers, one can derive both a practical and a theo- 
retical result: 

• An estimate of the quality of approximation one should have on each involved 
constant, such that the whole Boltzmann sampler is very unlikely to output any 
safety interval that does not contain the true value of the corresponding generat- 
ing function. Not surprisingly, @(logn) bits are enough when the expected size 
is n. 

• A (still hypothetical) construction for a truly exact Boltzmann sampler, if one 
assumes stronger oracles than those provided by [12]. If one assumes oracles 
that give both an upper an a lower bound for each involved constant, and that the 
oracles can be called repetitively to decrease the difference between these bounds 
(say, by a factor of 2 with every iteration), then the "safety interval" approximate 
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Boltzmann samplers can be modified into exact Boltzmann samplers which wiU 
call the oracles a small (logarithmic) number of times on average. 

7 NOT-QUITE BOLTZMANN SAMPLERS 

Boltzmann samplers share two important and useful properties for the user interested 
in practical random sampUng of large structures: (1) they output structures of random 
size, with the guarantee that any two structures of the same size have the same proba- 
bility of being obtained, and (2) they are "stable " under the constructions of Section 3. 

While property (1) is the essential one for direct practical applications, making it 
possible to use, say, rejection to transform Boltzmann samplers into exact- or approximate- 
size samplers, property (2) is the one responsible for the wide applicability of the 
method. Nevertheless, sometimes property (1) alone can be obtained while keeping 
algorithms of low complexity, typically by using a "nonstandard" final construction 
using Boltzmann samplers as subroutines. 

A good example of this is the Hadamard product of two combinatorial class. If 
^ and ^ are two arbitrary (unlabelled) combinatorial classes, their Hadamard prod- 
uct = ^ =^ is the subset of ^ x ^ that only contains pairs c = (a, b) with the 
same size (with |c| defined to be \a\ = \b\ instead of the sum for the classical product). 
The corresponding generating function is none other than the Hadamard product of 
generating functions, 

C{z)=AQB{z)=Y,anb„f, 

n 

with radius of convergence at least the product of radii of A(z) and B{z). 

A real Boltzmann sampler for can be written easily: one simply checks that, if 
X = xaxb with xa < Pa and xg < ps, the algorithm 
repeat 
a •<-r^(xA) 

p ^ Ym[xn) 

until |a| = |J3| 
Return (a,/3) 

terminates with probability 1 and in finite expected time, and outputs a structure 
under the Boltzmann distribution with parameter x. It is, however, inefficient: each 
iteration has success probabiUty C(x) / {A{xa)B{xb)^, which can be very small when x 
is close to the dominant singularity. 

Bodini et at. [3] use the classical Birthday paradox to devise a more efficient algo- 
rithm that preserves property (1), though not the whole Boltzmann distribution: simply 
alternate (either deterministically or randomly) between Fs/ (xa) and Y,'^{xb), retain- 
ing only the first obtained structure of each class and size, until a pair with the same 
size can be formed. 

8 CONCLUSION 

Boltzmann samplers are an attractive class of random generation algorithms for classes 
of combinatorial structures that lend themselves to combinatorial decompositions, al- 
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lowing for fast (linear-time, or quasi-linear-time) generation of structures with size in 
the millions for simple classes, and well into the tens of thousands for complex classes. 
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